Suggested references:
C Gardiner, Stochastic Methods: A Handbook for the Natural and Social Sciences, Springer.
D Higham, An algorithmic introduction to Stochastic Differential Equations, SIAM Review.
Statement of problem: given function $g \left( X^{(1)}, X^{(2)}, \dots, X^{(n)}, X, s \right) = 0$ defining an $n$-th order differential equation
$$ \begin{equation} \frac{dX}{ds} = f \left( X(s) \right) \end{equation} $$where we refer to $f$ as the drift function for reasons that will become clear later.
The formal solution to this equation is
$$ \begin{equation} X(t) = X_0 + \int_0^t f \left( X(s) \right) \, ds. \end{equation} $$Note a simple example. Given
$$ \begin{equation} \frac{dX}{ds} = a X(s) \quad \implies \quad X(s) = X_0 e^{a s} \rightarrow_{s \rightarrow + \infty} \begin{cases} \infty & a > 0 \\ 0 & a < 0. \end{cases} \end{equation} $$Solve this and think about how many solutions may exist.
Think about the integral form of the differential equation
$$ \begin{equation} X(t + h) - X(t) = \int_t^{t+h} f \left( X(s) \right) \, ds \approx \begin{cases} f \left( X(t) \right) \, h & \\ f \left( X(t+h) \right) \, h & \end{cases}. \end{equation} $$Picture here
This gives two numerical methods:
$$ \begin{align} X(t+h) - X(t) & \approx f \left( X(t) \right) \, h && \leftrightarrow & X_{n+1} - X_n & = f(X_n) \, h & \text{Explicit Euler}, \\ X(t+h) - X(t) & \approx f \left( X(t+h) \right) \, h && \leftrightarrow & X_{n+1} - X_n & = f(X_{n+1}) \, h & \text{Implicit Euler}. \end{align} $$Both these methods are first order in the sense that, given an approximate solution $X^h(t)$ and the exact solution $X(t)$, we have that
$$ \begin{equation} | X(T) - X^h(T) | \le C(T) h \end{equation} $$where $C(T)$ is a constant depending only on the time at which this is evaluated.
Apply the explicit Euler method to the simple problem
$$ \begin{equation} \frac{dX}{dt} = -a x, \qquad a \gg 1. \end{equation} $$This is a stiff problem.
The algorithm becomes
$$ \begin{align} && X_{n+1} - X_n & = - a X_n h \\ \implies && X_{n+1} & = X_n \left( 1 - a h \right). \end{align} $$This will only converge to the correct answer (0!) when $|1 - a h| < 1$ which implies $h < 2 / a$ (as $h > 0$). When $a \gg 1$ this is a severe restriction.
Compute the step size restriction for the implicit Euler method.
What are we trying to do? Consider the motivating example
$$ \begin{equation} \frac{dX}{ds} = \left( a + \text{"error"} \right) X(s), \quad X(0) = X_0. \end{equation} $$The "error" should average to zero.
More formal definition is
$$ \begin{equation} X(t, \omega) = X_0 + \int_0^t f \left( X(s, \omega) \right) \, ds + \underbrace{\sigma B \left( t, \omega \right)}_{= \sigma \int_0^t dB(\omega, s)}. \end{equation} $$Each time I repeat the experiment I get a different answer. The first part is "deterministic" and looks like the ODE; the latter part is new.
Now, we list some assumptions of varying strengths.
Assumption 1. It should be clear that $B(0, \omega) = 0$ (no "error" accumulation at the initial time).
We can also look at two close points in time to see that
$$ \begin{equation} X(t+h, \omega) - X(t, \omega) = \int_t^{t+h} f\left( X(s, \omega) \right) \, ds + \sigma \left( B(t+h, \omega) - B(t, \omega) \right). \end{equation} $$Assumption 2. Given our discussion we assume
$$ \mathbb{E} \left( B(t+h, \omega) - B(t, \omega) \right) = 0. $$We now drop $\omega$ from the notation to make explicit
Assumption 3. An independence assumption:
$$ \begin{equation} B(t_2) - B(t_1) \underbrace{\perp}_{\text{independent}} B(s_2) - B(s_1) \quad \text{for all times } s_1 < s_2 < t_1 < t_2. \end{equation} $$Assumption 4. Final independence assumption: $B(t+h) - B(t)$ depends only on $|h|$ and not on $t$.
The last two assumptions (assumptions 3 and 4) suggest that this is a "pure" error term - that it has no predictable structure.
The error term, or Brownian motion, or Wiener process we will define satisfies:
We write this as
$$ \begin{equation} \int_0^t G(s) \, dB(s). \end{equation} $$This is slightly odd as we're integrating against a random process. What we're aiming for is a way of solving (in the integral form) the SDE
$$ \begin{equation} \frac{dX}{ds} = f(X(s)) + g(X(s)) \frac{dB}{ds} \end{equation} $$Split the interval $[0, T]$ into segments length $h$. On each segment we want the increment to be a normally distributed random variable, $N(0, h)$.
In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
In [2]:
T = 1.0
h = 0.01
N = int(T/h)
t = np.linspace(0.0, T, N)
dW = np.sqrt(h) * np.random.randn(N)
W = np.cumsum(dW) - dW[0]
In [4]:
plt.plot(t, W)
plt.xlabel(r"$t$")
plt.ylabel(r"$B(t)$");
Out[4]:
Rerun this to see that it changes. Rerun again with a smaller $h$. You'll see it gets "spikier". In the limit as $h \to 0$ the function is not differentiable.
In [6]:
M = 5
dW = np.sqrt(h) * np.random.randn(N, M)
W = np.cumsum(dW, axis = 0) - dW[0, :]
plt.plot(t, W)
plt.xlabel(r"$t$")
plt.ylabel(r"$B(t)$");
Try to solve
$$ \begin{equation} \int_{0}^T B(s) \, dB(s) \ne \frac{1}{2} \left( B^2(T) - B^2(0) \right) \end{equation} $$which would be the "normal" calculus result.
Instead do the "natural" (or naive) partial sum thing:
$$ \begin{equation} S_n = \sum_{i=1}^n B(\tau_i) \left( B(t_i) - B(t_{i-1} \right), \qquad t_{i-1} \le \tau_i \le t_i \end{equation} $$This is a random variable, as $B$ is a random variable. So this has a mean, variance, higher moments...
This means that, depending on how you choose your $\tau_i$, you will get different answers for the exact solution!
where $a$ is the relative distance between $t_{i-1}$ and $t_i$ at which $\tau_i$ is located, using the earlier property $< B(t), B(s) > = \min\{t, s\}$.
The Ito integral uses $a = 0$ so that $\tau_i = t_{i-1}$. This implies $< S_n > = 0$.
We are trying to solve
$$ \begin{equation} X(t) = X_0 + \int_0^t f(X(s)) \, ds + \int_0^t g(X(s)) \, dB(s). \end{equation} $$The simplest thing to do is to use the Ito integral applied to the final term. By our independence assumptions, as the expectation of the Brownian path is zero, the expectation of the final integral (our "error term") is zero, as we want.
So the equivalent of the explicit Euler method is
$$ \begin{equation} X(t+h) - X(t) = h f(X(t)) + g(X(t)) ( B(t+h) - B(t) ). \end{equation} $$The increment in the last part is random, as it's a random variable!